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I. INTRODUCTION 



While the values of the critical parameters for the fully symmetric hard-sphere model 
of an electrolyte (the "restricted primitive model" or RPM have received extensive 

theoretical attention for a number of years (see, e.g., P-|lO|]), only recently has interest 
focussed on the effects of asymmetry on phase-coexistence and criticality [[TTHT^. Sabir, 
Bhuiyan, and Outhwaite used the mean spherical approximation (MSA) and two dif- 
ferent Poisson-Boltzmann approaches to compute critical parameters resulting from both 
size and charge asymmetry; the reduced critical density, p*, was always reported to increase 



were 



14 and 



with greater size asymmetry, but the trends for the reduced critical temperature, T*, 
inconsistent. (For appropriate definitions of reduced temperature and density, see 
below.) Gonzalez- Tovar used the MSA and found, via the energy route, that both T* and 
p* increased with greater size asymmetry ||12| . Most recently, Raineri, Routh, and Stell also 



employed the MSA, but augmented the analysis by incorporating Bjerrum-Ebeling-Grigo 
pairing; they likewise predicted that both T* and p* increase monotonically with size asym- 
metry HTSf. Recent simulations by Romero- Enrique et al. |]14| (see also [l^), however, reveal 
the opposite trends — namely, that critical temperature and density decrease strongly with 
increasing size asymmetry. Hence, as remarked in [0], our current state of even qualitative 
theoretical understanding appears less than adequate. 

The present report therefore formulates and analyzes extensions of the Debye-Hiickel 
(DH) approach to asymmetric primitive models [|T7|. We may recall that the original DH 



theory, when supplemented by Bjerrum ion pairing and dipolar-pair solvation in the ionic 
fluid 1^, remains the most quantitatively successful theory of coexistence and criticality in 
the size-symmetric RPM (see, e.g., [7(b), 10]). One might, thus, reasonably hope that a 
suitable extension of the DH and (DHBjDI [Q) theories to unequal ion sizes will, at least, 
yield correct trends in the dependence of T* and p* on asymmetry. 

In fact, in their original 1923 paper, Debye and Hiickel [|I|] already claimed to have explicit 
results for asymmetric primitive models (see ^\)- However, serious flaws in this historic 
formulation suggest other, more systematic extensions of symmetric DH theory appropriate 
for the asymmetric case [|T7|. We analyze these improved "asymmetric Debye-Hiickel" (ADH) 
theories by comparison with exact series expansions and newly developed bounds, and also 
assess their general physical character. We find, in fact, that, contrary to Refs. |TT|-|T3[, the 



critical parameters predicted by our modified ADH theories [0 exhibit trends in agreement 
with those obtained in the recent simulations 
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In order to focus on the effects of size asymmetry, we confine the present discussion 
primarily to the two-species size-asymmetric primitive model (or SAPM), which consists of 

equal numbers, A^_|_ = A^_, of positive and negative ions with hard-core diameters a > a++, 

and charges of equal magnitudes q+ = — g_ = z±qo. (Of course, the complementary case 

a++ > a follows trivially by symmetry.) All material and space is assumed to possess a 

uniform dielectric constant, D. For the most part, we also assume additivity of diameters, 

a_|__ = a = |(a++ -|- a ). (1.1) 

The degree of asymmetry will be described either by the fractional deviations from a, namely, 
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or by the diameter ratio WM, 



A = a__/a++ = (l + (5_)/(l-(5_), (1.3) 

where we have used the fact that 5+ = — 5_ for the additive case. Note that the general 
asymmetric primitive model electrolyte possesses a large — in fact, three-dimensional — 
parameter space: two degrees of freedom for size asymmetry (since additivity need not be 
assumed) and one for charge asymmetry (noting that the overall magnitude is irrelevant). 

We take the total number density to be p = N/V, with V the volume and = + A^_. 
The reduced density is then appropriately defined by |^ 



p* = {p+ + p.)a\ (1.4) 

which seems the most relevant parameter since the critical densities are typically low, p* ^ 
0.1. But we note that other workers have used alternative definitions |1TT|JT3[| . The natural 
temperature scale for studying criticality (i.e., low temperatures) in the SAPM is set by the 
energy of an attractive pair of ions at contact, so we take 

T* = kBT/i\q+q_\/Da^_) . (1.5) 

It is also convenient to define the reduced electrostatic energy, normalized by the energy of 
closest approach of a (+, — ) pair, namely, 

u{p,T) ^ U'^ip,T)/Ni\q+q^\/Da), (1.6) 

where (= Un — ^NksT) is the overall excess energy. 

The theoretical description of size-asymmetric models is, naturally, more complicated 
than in the symmetric case. Specifically, in considering the correlations and fluctuations in 
the neighborhood of a particular ion (which, following DH, we may suppose is fixed at the 
origin), three distinct surrounding spherical shells or zones must be accounted for. To see 
this, suppose first that the selected, fixed ion is positive and that the diameters are ordered 

in size according to a+_|_ < a_|__ < a : we term this the "inner" case since the like-like 

diameter of the central charge is smaller than a+_. Trivially, no ion center can enter the 
"interior" zone < r < a++, where r is the radial distance from the origin; of course, this 
applies to the symmetric RPM as well. New to the SAPM, however, is the "border" zone 

(or, more precisely, the "stx^-border" zone) a++ < r < for this "inner" scenario, that 

is shown shaded in Fig. |l|: this can be populated only by the centers of the smaller, i.e., 
positive, ions; the larger negative ions are excluded. Finally, as in the symmetric model, 
positive and negative ions may be present in the exterior zone r > a+_. When a larger 
negative ion is chosen to be at the origin, there is a complementary "super-border" zone, 
a+_ < r < a , into which only positive charges can penetrate. 

The presence of these "charge-unbalanced" zones in the asymmetric case turns out to play 
a crucial role which, in particular, means that the simple extension of the DH treatment 
(ADH theory), although in agreement in leading orders at high T with the exact series 
expansions P,p!6|JT7[], fails badly at lower temperatures. This is shown below in Sec. 11 after 



we develop the ADH formulation, which allows for the border zones. While the failure 
is not entirely surprising, since the general DH formulation relies on a high-temperature 
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expansion, it does present a sharp contrast to the symmetric case. One striking consequence 
of the difficulties is the theory's violation, at low temperatures, of an extension of Onsager's 
lower bound for the (internal) energy of the RPM that we establish in Appendix A for 
general primitive models that are asymmetric in both charge and in size. 

It is an interesting fact that analysis of the general primitive model with nonadditive 
diameters presents few additional problems in most cases. (See, e.g., and references 
therein.) On the other hand, the nonadditive models are potentially useful in applications 
because they can represent short range interactions beyond the pure Coulombic couplings. 

Thus if, for example, a++ = a = oq, so that one has a size- symmetric model, but a+_ > oq, 

a moment's consideration shows that unlike ions experience a strong short-range repulsion 
relative to the geometric ion size oq. This competes with the ionic attractions and, thus, for 
example, increases the relative size of a Bjerrum pair and decreases its stability. Conversely, 

the case < oq can be viewed as representing an enhanced short range repulsion between 

like ions of a "true" diameter less than oq. In reality, of course, the interactions in ionic 
systems deviate from pure electrostatics at short distances. Accordingly, although we mainly 
focus below on the additive case, we briefly indicate where necessary how to construct DH 
theories for the general, nonadditive model [1^. We believe this theory is of interest in its 
own right as well as providing a basis for more elaborate treatments. Further details of the 
linear asymmetric DH theory are also given in ||l7 . 



Because of the failures of the linear asymmetric DH theory, we have explored some 
modifications that, for modest degrees of asymmetry, ensure satisfaction of the energy bound 
of Appendix A and of thermal convexity requirements [0; these criteria are also satisfied 
for large asymmetry in the critical and coexistence regions. The modified ADH theories 
accomplish this by specifically allowing for, and compensating in a natural way, the charge 
imbalances induced by the existence of the border zones [|r^. These theories are expounded 
in Sec. HI and their predictions for criticality are examined in Sec. IV. It transpires, as 
mentioned above, that the simplest modified theories appear to be the ffist theoretical 
treatments to predict correctly (as judged by the simulations [0,|1^) that both T* and p* 

fall when size asymmetry is introduced; however, for A = a /a++ ^ 4, the critical density 

displays incorrect, nonmonotonic behavior. Finally, Sec. V presents a brief overview and 
some conclusions. 
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II. ASYMMETRIC DEBYE-HUCKEL THEORY 



A. Formulation of the theory 

In this subsection we present the ADH theory. Since the cost in additional complexity 
is slight, we consider here (and in Appendix A) a general two-species model with point 
charges = z^q^ (a = +, — ) centered in hard spheres of diameters a^^^ with the collision 

diameters, a„r, restricted only by a++ < a = a+_ < a The Debye-Hiickel procedure [|I|, 

interpreted generally, directs one to calculate the excess thermodynamic functions, based 
on approximations for the correlation functions which are then integrated via the "energy 
route" . Specifically, defining ip^ to be the average electrostatic potential at a (fixed) ion 



of species a due to all other ions, the reduced electrostatic energy, defined in (|T]^), is given 

by 

where the superscripts < and > merely serve as reminders of the relative ion sizes, and also 
facilitate formulation of the nonadditive case. Other thermodynamic quantities, such as the 
free energy and pressure, follow in principle after suitable integrations and differentiations. 
We derive an explicit closed form expression for the internal energy as a function of T and 
p, and the second virial term for the pressure is also discussed. 

We begin the calculation of the ion potentials ipaiT,p) by fixing an ion of species a at 
the origin, as in Fig. ^ The induced electrostatic potential, 0(r), and corresponding charge 
density, Pq{r), are then related by the (exact) averaged Poisson equation, namely, 

V^(0(r)). = --(p,(r))., (2.2) 

where the subscript a indicates that the averages are performed with a charge of species a 
at the origin |]l9l . The ion potentials follow from the limit 



tpa{T,p) = lim 



Dr. 



(2.3) 



which eliminates the self-interaction of the fixed charge at the origin. 

The DH procedure then introduces two approximations for the average charge density, 
{pq{r))cr- The well-known Poisson-Boltzmann approximation is followed by a linearization 
of the exponential in the Boltzmann factors, yielding 

{Pq{r))a^ (lrPrexp[-pqr{(f){r))^] , (2.4) 

T= + ,- 

~ 53 g.p.[l-/5g,(0(r)).] , (2.5) 

T=+ ,- 

where (3 = l/ksT. In the size-asymmetric models, the approximate charge density must be 
allowed to take a different form in each of the three distinct zones around the central charge. 
Thus, when the fixed charge is positive {a = +) one has 
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{pqij))+ = (l+^ir) for r < a++ , (2.6) 

- [1 - for a++ < r < a , (2.7) 

~ -(Ky47r)((/)(r))+ forr>a, (2.8) 



where, as usual |T9|, overall electroneutrality has been imposed. Here, the inverse Debye 



length is defined in the standard way via 

{Koaf = x\ = {Ana^ql/DkBT) ^ p^zl . (2.9) 

Note: (i) there is no approximation in the innermost exclusion zone (r < a++); (ii) the 

sub-border zone (a++ < r < a = ) can contain none of the larger negative ions so that in 

( p.7[ ), which represents the fundamental extension of the original DH theory, only g+ appears 
on the right-hand side; and (iii) the exterior zone (r > a) follows the standard Debye-Hiickel 
form as in the symmetric RPM |ll.T9[ . 



The approximations (|2.6|) - (|2.8|) complete the reduction of the averaged Poisson equation 
( p.2| ) to the basic ADH equation for a + ion. This may be solved for four unknown coefficients 
using standard electrostatic boundary conditions (continuity of and d(f)/dr) at r = a++ 
and r = a_|__ = a with 0(r) ^ as r — oo |]19|, from which ip^ follows via (|2.3|) . The closed 



forms for the electrostatic energy in the ADH approximation then result from combining the 
general expression for the energy ( p.lD with the appropriate expressions for the ipa potentials, 
which are presented below. If we introduce 

9,^[\z^\/{z+-z^)f\ (2.10) 

one finds that the potential in the "inner" scenario, with a smaller central positive ion 
(indicated by the superscript <), may be written 

V'+ (go; T, p- {z^l, {(5^1) = j — , (2.11) 

where the various contributions are 

A< = {1 + xd + xdS+) [cosh {e+XDS+) - 1] , (2.12) 

^2 = -iK^ + + d+XDS+) sinh {e+XD6+), (2.13) 

A< = xd6+, (2.14) 

B"^ = xo cosh {9^xr)S+) — O+xo sinh {9+xdS+), (2.15) 

C< = {I + xd + xdS+) cosh {e+XDS+) - 9+{e^^ + XD + xd5+) sinh {e+XD5+). (2.16) 

Note that 5+ < for this < case. The A^ terms of ip'^ have been grouped so that each 
vanishes individually when 5+ 0; in this limit of size symmetry, furthermore, the potential 
reduces to the standard DH result |I],p|,p!9| for a size- symmetric primitive model, namely, 
V'^ = V'+ = —{z+qQl Da)xj:,/ {1 + X£)). When the diameter of the negative ion is less than a, 
one should simply switch the species subscripts in ( |2.11| ). 

For the "outer" situation, with a larger central ion (indicated by superscript >), the 
corresponding result is 
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^lj>[qo]T,p-{z,},{5„}) = I 1 ^ 



where the contributions are now 

A> = {1 + xd) [cosh {9+xd5-) - 1] , (2.18) 

A> = {0-^ + O+xd) sinh {e+XD6^), (2.19) 

A> = -xd5-, (2.20) 

= Xd cosh {9^xdS~) + O^xd sinh (Oj^^xdS-), (2.21) 

C> = (1 + Xd) cosh {e+XD6-) + {6^^ + O+xd) sinh (O+xdS-). (2.22) 



Naturally, ?/^^ also reduces to the symmetric DH result when 6-^0. (If a++ > a , one 

should simply switch the + and — subscripts.) 

For future reference, notice that both ip^ and ip^ have readily-calculated zero- 
temperature limits, namely, 

,<(^.o) ^ - (^) ^ •^?(^-°)-^. 

which are independent of density. Furthermore, small diameters, a, are of interest; but ip^ 
diverges for "point" ions, i.e., in the limit ag-o- —>■ (or 6^^ —>■ —1). This, in fact, serves as a 
warning sign, as will be seen below. 



B. Assessment of ADH Theory 

We now examine the predictions of the linear ADH theory. Following Debye and Hiickel 
P|JTP|], the free energy is to be obtained from the internal energy in ( |2.1D by employing the 
Debye charging process (which is equivalent to using the standard thermodynamic relation 
for the free energy in terms of the energy). This yields the reduced excess free energy density 
as 



r{T, p- {z^}, {5J) = -A^^{V- T)/VkBT 



dC ^+(Cgo)-^>(Cgo) 



(2.24) 



where is the total excess Helmholtz free energy. It must be recalled that in addition 
to the explicit dependence on qq (and, hence, on () entering (^]l]) via (|2.11| ) and (|2.17|) , an 
implicit dependence occurs via kd oc go which enters (|2.12| )- (|2TT^ ) and (|2.18|) -( pr22D . 

Except in the standard symmetric case (A = 0, 6+ = 6- = 0) most of the integrals 
involved in ( |2.24| ) seem intractable. Nevertheless, one may derive an expansion for the 
free energy, and thence for all other thermodynamic quantities, such as the pressure, by 
expanding the expressions ( |2.11| )-( |2.22| ) and integrating term by term. By this route, we can 
check the theory against exactly known expansions (see, e.g., |^,16,17|). 

Thus, consider the low-density/high-temperature expansion for the general primitive 
model (with arbitrary charges q„ = z^q^ and diameters a^r) which is known to overall order 
p^/^ in the density. On specializing to the two-component case, this may be written 
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CO 

+ E 52,n/T*" + 0[p'/' In (/tz^a)] , (2.25) 

n=0 

where we again recall Kjya = xd oc a/p*. The leading term here varies as p^/^ and represents 
the DH limiting law reproduced by all sensible approximations. The coefficient i?2,o derives 
from the second virial coefficient for a pure hard-sphere gas. At the present level of approx- 
imation, this may be accounted for precisely by including the hard-core free energy density, 
/^*-^, in the basic approximation The remaining coefficients -82,1, -82,2, -B2 3, -82,3, -82,4, • • • 
arise from the electrostatic interactions. 

Now, DH-based approximations (in common with the MSA, etc.) cannot generate the 
p^ (In p)/T^ term with coefficient -82,3- However, our linear ADH theory yields the exact 
leading coefficients in order 1/T and Explicitly we find 



B21 = 7ra^ \ z^z 



■5+(2 + 5+) + 5_(2 + 5_ 



[z+ - z_Y 

which vanishes for the size-symmetric models (with (5o- = 0), and 

zibs ^zH/ 



(2.26) 



zV 



(2.27) 



It is, indeed, striking that the exact leading order dependences in powers of 1/T are re- 
produced for arbitrary diameters, a++, a+_, and a (including nonadditivity and point 

charges). However, one probably should not be so surprised since the arguments leading 
to the basic ADH equations, ( |2.2| ) with ( ^.4[ ) and ( P3| ) are, on reflection, clearly valid to 
leading orders in p and 1/T even if it is not obvious how far the validity will go. 

To study the ADH theory more quantitatively, we focus on the energy for the additive 
1:1 case (2+ = — z_ = 1). No integrations are then required. Fig. ^ displays energy isochores 

as functions of T* at high and low temperatures for degrees of asymmetry A = a /a++ = 1 

(the RPM for reference, solid curves), 2, 4 and 6. These values of A correspond to size 
deviations = 5_ = 0, |, | and |, respectively. It is evident that ADH theory predicts that 
size-asymmetry lowers the electrostatic energy, with the effects being greatest at the highest 
and lowest temperatures. At intermediate T the predictions for varying A are grouped 
together according to density. 

Notable features of the high-T plots of Fig. ^(a) are, ffist, that the finite values of u 
attained when T ^ 00 are exact. This can be verified by a simple a priori calculation; 
alternatively, the values may be checked from the exact order p^/T term as displayed in 
( p.25| ) and ( p.26|) that is correctly reproduced by ADH theory. The limiting slopes at T* = 00 



are discussed in further in |T] 

The behavior of u{T,p) at low temperatures proves, however, much less satisfactory. 
From Fig. ^(b) one sees that the energy is predicted to fall increasingly rapidly for increasing 
asymmetry when T falls, passing well below the ground-state prediction u = — ^ of the 
symmetric DH theory; but see also for the DHBjDI extensions [^]. Most seriously. 



however, for large enough asymmetry, the energy drops below the 1:1 size-asymmetric lower 
bound, namely, u > —1, established in Appendix A (which also treats general z^ and z_). 
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More concretely, the ground states predicted by ADH theory follow from ( p.l|) and (p.23|) 
which yield 



,ADH 



u 



(T = 0;A) = -i(3 + A). (2.28) 



This expression evidently violates the bound when A > 5 ( or = 5_ > |). 

Of course, this is unacceptable, but even for smaller values of A this behavior casts doubt 
on the value of the approximation at low T in the vicinity of the expected critical region. 
The origin of this behavior can be understood physically by examining the ADH predictions 

for the mean total charge, say Q^, within the sub-border zone a++ < r < = a which 

can only be penetrated by like, i.e., positive charges. Recalling the notations in ( |2.2| ) and 
the approximations (|2.4| )- (|2.8| ), the theory yields 



a 4 



Q< = / d'r (p,(r)) + 

Ja++ 

- P"dVg+p+ [1 - /?g+(0(r))+] , (2.29) 
where the precise form of (0(r))_|_ follows from the explicit solution of ( |2.2| ) with ( p.6| ) 



17| . Now, when T oo the approximation becomes exact, yielding Qf = g+p+V^, where 



= |7r(a+_ — ct++ ) is the volume of the sub-border zone (see Fig. 1). At finite but 
large temperatures, ( p. 291) correctly predicts the linear decrease of with P as the (+, +) 
repulsions come into play. Of course, the repulsions become increasingly effective as T 
decreases so that, at low T, one expects to remain positive but to be of magnitude only 
F^g+p+ exp (— /?g^/i5a+_) (although, because, of the negative screening cloud that builds 
up for r > a+_, this estimate might need to be modified somewhat). 

However, if one neglects the screening, the mean potential (0(r))_|_ in (|2.29| ) will be of 



order +q^/Da^^, or larger, in the sub-border zone. Then, as T falls and /5 ^ oo it is 
evident that ( |2.29|) is likely to predict, first, a totally unphysical negative value for (even 



though no negative charges can enter the sub-border zone) and, eventually, a divergence of 
to — oo. Explicit calculations fully bear this out. For example, when A = 2, and 



p* = 0.01, one finds negative for T* < 0.9; for p* = 0.1 (> p*) the change of sign is 
delayed until T* ~ 0.6 but the divergence to — oo is more rapid. Since the ADH theory 
for the border zones embodied in ( p.5| ), ( |2.7| ), and ( |2.29| ) reflects a high-T expansion, the 
problems at low T should not be a great surprise. On the other hand, the successes of the 
original DH theory for the symmetric case at low T (e.g., P,p!o[|) clearly hinge on the absence 
of any "charge-unbalanced" border zones. 

Clearly the serious physical defects of the ADH treatment must be rectified if the theory 
is to have any value for T* ^ 1: and, we expect, T* ^0.1 ||7|-[T3[]. Before addressing that 
task, however, we present, briefly, the original claims of Debye and Hiickel fl]] as regards a 
system of ions with arbitrary diameters. 



C. Original DH Theory 

In their original 1923 paper, Debye and Hiickel did, in fact, discuss the case of an 
"arbitrary ionic solution" — that is, within the model they adopted, a mixture of hard- 
spheres of varying diameters Oo-o- and charges q^, or the general primitive model. However, 
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in their analysis no mention is made of the border zones in which charges of only species 
with small enough diameters can be present. (Note that in the multispecies case there 
will, in general, be a number of distinct border zones around each charge.) Rather, their 
argument is presented as if all species were of the same diameter |lT|,p!7[|. The mathematics 



is then identical to that for the RPM. In our notation, Debye and Hiickel thus present the 
conclusion 

= - i¥) TT^' P.30) 

\ U / 1 + Undaa 

(i.e., precisely the RPM result when a^^^ = a). The equivalent free energy for a 1:1 model is 
then 



IGtt [ (1 + 5+)3 

, 2 In [1 + (1 + 5^)xn] + [(1 + 5^)xd]^ - 2(1 + 6^)xd \ 



from which predictions for criticality etc. follow: see below. 

Neglecting all the border zones spares this "original DH" theory the pathology of our 
ADH theory, but seriously reduces the plausibility of their treatment for size-asymmetric 
models. Thus, not so surprisingly, the high-T/low-p expansions of the thermodynamics 
resulting from ( |2.31 ) differ at the first correction to the limiting behavior — see i?2,i in 



( ^.2(j| ) — from the exact results captured by our ADH theory; indeed, the former depends 
only on the ratio x'jj oc p/T, while the exact results are more complex. 
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III. MODIFIED ASYMMETRIC DH THEORY 



We now discuss a class of modifications of ADH theory (originally introduced in |r7[ ) 
which avoid the unphysical behavior of the charge density in the border zones: this be- 
havior was identified in Sec. II. B above as a root cause of the pathological low-T behavior 
which included violations of the lower bound on the internal energy (Appendix A). In the 
following section, we examine the predictions of these modified ADH theories for the critical 
parameters (including a comparison with the original, 1923 DH proposal for asymmetric 
systems) . 



A. Introduction of Border Zone Factors 



It is clear from the previous considerations of the ADH prediction ( p.29|) for the charge 
in a sub-border zone, that the standard DH linearization of the Poisson-Boltzmann expo- 
nential is generally inadequate in any zone where the mean charge is necessarily unbalanced 
because of strongly repulsive short-range, steric repulsions. As an alternative first step, that 
will still yield an analytically tractable theory, we forego the self-consistent aspect of the 
Poisson-Boltzmann approximation ( p.4|) in the border zones and consider replacing the self- 
consistently determined electrostatic potential {(f){r))„ by a fixed, but possibly temperature- 
and density-dependent effective potential, (p'^r^T, p). 

One possibility is merely to use the bare potential, i.e., to put (pl = Qa/Dr. The right- 
hand side of ( |2.7] ) would then read exp (— /3g^/Dr); but an obvious drawback of such 
an approximation is that the resulting closure of Poisson's equation, (p.2|) , is no longer 
readily solvable. Instead, we simplify further by neglecting the r-dependence of Thus 
we introduce temperature-dependent border zone factors J-''^{T) and J^^{T) and consider 
the replacement of (|2.7|) for a sub-border zone by 



{p,{r))+=p<{r-T,p)^q+p+T<{T), 



< r < a 



(3.1) 



Likewise, for a super-border zone (indicated, as before, by the superscript >) we advance 



a < r < a 



(3.2) 



The calculation of the tpa potentials now proceeds just as in Sec. || but the results take 
a simpler form, which we distinguish by using a circumflex. Maintaining our convention of 
a smaller positive ion (5+ < 0) and larger negative ion (5_ > 0), we find 



1 + Xd^t: ^ 

Z+T*6-J^>{T) r 



Da \ + xd 
Da 1 + xd 



2x2 ■ 



2 + 5+- xd{5+ + U_ 



1 + Xd 



2iz, 



2 + 5_+xd{5- + \61) 



(3.3) 
(3.4) 



These modified potentials, valid for general J-'{T), reproduce the symmetric DH theory 
when 5o- — > 0, and always generate the standard limiting laws when p — > 0. For large 
xn oc (p/T)^/^, these approximations for the -ipa- behave as x'jjTJ^, whereas the standard DH 
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expressions approach constants. While the consequences of this fact are not obvious, some 
numerical results are discussed in Sec. III.C, below. 

On the other hand, both of the high-T second- virial terms, -82,1 and -82,2 that are correctly 
generated by the ADH theory [see (p.26|) and (|2.27|) ], will be reproduced now if the zone 
factors satisfy 



1 - 

1 + 



\z+/z- 



T*(l + i5+) 
1 

T*(l + i5_) 



+ 
+ 



1 



1 



(3.5) 
(3.6) 



In principle, JF< and JF> could be chosen to generate further exact second-virial terms; but 
such an approach does not seem of much practical value. 



B. Choice of Zone Factors 

Undoubtedly the simplest reasonable approximation for the zone factors is provided by 
taking 

MF: J^<(r) = T>{T) = 1, (3.7) 



or, equivalently, by dropping the term linear in j3 in ( |2.5| ) and ( |2.7] ). The results will fail 



to reproduce the correct -82,2 coefficient in (|2.27|) , although the exact i?2,i term will be 



generated. Physically, this approximation amounts to a direct mean-field or van-der-Waals- 
type approach in which the effects on the internal energy of the (T = 00) unbalanced zone 
charges, = q+p+V^ and = q^p^V^ (where and are the zone volumes), are 
accounted for in a direct way that neglects fluctuations. The zone charges remain flnite and, 
in fact, fixed for all T: but, at least for small relative zone volumes V^/a^ ~ 47r(5+(l -|- 6^) 
and /a^ ~ 47r5_(l + 5_), one might reasonably expect that the initial thermodynamic 

trends with increasing asymmetry A = a for small A, will be correctly predicted. 

A second natural step, following our initial discussion, is to take ^''"(r) to be the direct 
potential, q„/Dr, evaluated at the mean radius of the border zones, namely, r< = |(a++ -|- 
a) = [1 + |(5+)a and = {1 + \5-)a. This amounts to adopting 

EXP<: ^<(T) = exp [-pql/Da{l + \5+)] , (3.8) 

and similarly for J^^{T) but with the exponent 

y> = +l3q+q^/Da{l + \5.) . (3.9) 

It is encouraging that these choices precisely satisfy ( |3.5| ) and ( p.6| ) and so reproduce both 
S2,i and i?2,2- 

As regards the sub-border zone where, as discussed in Sec. II. B, one expects (5+(T, p) to 
vanish at low T (or, at least, approach very small values), the assignment (|3.8| ) seems rather 
satisfactory. Indeed, combining ( |3.1D and the flrst part of ( |2.29D shows that will remain 
positive but decrease exponentially rapidly as the (+, +) repulsions "flush out" charge from 
the sub-border zone. 
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It is worth pointing out that the exponential treatment of the sub-border zone can be 
extended within the DH self-consistent spirit by linearizing the Poisson-Boltzmann factor in 
( p.5| ) about the central value of the direct potential in the sub-border zone. Thus, if —y"^ 
denotes the exponent in (|3l8| ), while y = /?g+(0(r))_|., one may replace (|2.7|) by 

^q+p+e-y''[l-{y-y<)] . (3.10) 

The overall factor e~^^ should now ensure the sensible behavior of the sub-border zone 
charge, Q^, while the linear factor y oc {pg{r))^ accounts self-consistently for some fluc- 
tuation effects while still allowing integration of the closed ADH equations. We have not, 
however, examined this approach further. 

In the super-border zone case, however, matters are significantly different. Certainly, as 
follows from ( |3.6| ), one expects the mean border zone charge Q^{T, p), which will be positive 
for T = oo, to increase initially when T falls. This is simply because the Boltzmann factor 
enhances the attractions between the larger, central negative ion and the smaller positive 
ions. Nevertheless, the growth of cannot continue indefinitely as would be implied by 
adopting the EXP^ form with ( |3.9| ) for JF>. Regardless of other effects, the hard cores of 
the -|- ions must limit the number that can be packed into a super-border zone at any T: 
hence must saturate. This effect may be incorporated into the present framework via a 
saturation fraction, JF^ = s, which prescribes the limiting T — > super-border zone charge 
density as sq+p+] recall ( |3.2|) . A plausible zone factor is then the "regulated" exponential 
form 

REGEXP>: J^>(T) = — , (3.11) 

with c = s/{s — 1), which satisfies ( WM and approaches s for large (i.e., low T). In the 
limit s — > 1 (no growth in Q^) this reduces simply to the mean- field form ( |3.7D . 

The difficulty in utilizing this proposal, however, is to know what saturation value (s > 0) 
is appropriate. At low T, considerations of ion association indicate that the predicted 
saturation charge in the super-zone, namely, 

Q>, = V>q+p+s = Airqo "\ ^, '^'^ -' z+p*s , (3.12) 

1 + \z+/z_\ 

should not exceed the neutralizing value = l-Z-l^o- To meet this condition requires s oc 
1/p. While a density-dependent s and, hence, can be accommodated without changing 
( ^.4[ ), the subsequent expressions for the free energy become more complex. Furthermore, 
when p* increases, s must decrease and may even fall below unity. In these circumstances 
and in the absence of other effective selection criteria we believe it is appropriate to accept 
the mean-field form (p.7|) , i.e., to set JF> = 1 even though 52,2 will not then be correctly 
reproduced. 

In passing we mention, nonetheless, that we have also explored algebraic forms which 
respect ( p.6|) ||20|| . As an example, one can consider 



J^>{T) = {l+y>)/[l + {y>/s)r , (3.13) 
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which saturates at JF> = s (not necessarily greater than unity) when = | but decays to 
at low T when i/ > |. For s ^ 2 such approximations have an added, unexpected feature, 
namely, they tend to predict an additional asymmetric critical point at higher densities that 
does not smoothly connect to the standard DH critical point when A — > 0. While, owing 
to the strong short-range interactions implicit in the asymmetric systems, such extra critical 
points are not obviously unphysical (especially when the diameters are nonadditive), we will 
not discuss the algebraic forms ( |3.13| ) (or EXP^) further here. 



C. Behavior of Energy Isochores 

Before examining the predictions of the modified ADH theories for criticality, it is in- 
structive to examine energy isochores, such as those plotted in Fig. 3, for several densities 
near the critical value: compare with the results for the simple ADH theory in Fig. 2(b) but 
note the difference in vertical scales. At temperatures in the coexistence region (T* < 0.2) 
and for modest degrees of asymmetry (A ^ 3 - 4), both the MF-ADH theory [specified by 
(^] and the EXP<MF>ADH theory [where has been replaced by (^] predict that the 
internal energy increases above the RPM value with increasing asymmetry. This reduction 
in thermodynamic stability suggests, as we will confirm, that the predicted values of T* fall 
as A increases. 

We also find that the MF-ADH and EXP<MF>ADH theories do, indeed, satisfy the 
energy bounds of Appendix A for all (p, T) states relevant to the critical and coexistence 
regions in the additive case. Bound violations can occur, but these arise only at the highest 
densities (p* ~ 1) and when the asymmetry is great (A ^ 1). While this behavior under- 
mines the two modified ADH theories investigated as overall descriptions of the asymmetric 
primitive model, the pathologies occur far from the critical region and at unphysically large 
asymmetry levels. 

For moderate to large asymmetry and high density, however, we find that the MF-ADH 
and EXP'^MF^ADH isochores exhibit nonmonotonic behavior in T, indicating a thermo- 



dynamic instability as previously found in a variety of ion-pairing theories ||10[- For both 
the MF-ADH and EXP^MF-^ADH theories, these convexity violations occur for additive 
models (5_ = —5+), roughly, only when 6- + p* ^ 1. No convexity violations are found for 
A ^ 2 at any physical density (recalling the packing limit): this confirms the view that the 
present theories are of greatest validity for modest asymmetries. 
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IV. PREDICTIONS FOR CRITICAL PARAMETERS 



The discussions presented above of the various proposed modifications of the asymmetric 
DH theory indicate that simpler versions should provide a reasonable basis for further ex- 
ploration, at least in the case of small size asymmetries. Accordingly, we present here, first, 
by way of calibration, the predictions of (i) the original, 1923 DH theory embodied in ( p. 311) ; 

(ii) the MF-ADH theory which uses the simple, mean-field border-zone factors ( p. 71) ; and 

(iii) the EXP^MF^ADH theory which retains the mean-field treatment for the super-border 
zones but recognizes, via (|3.8|) , the decrease in sub-border zone charge resulting from the 
Boltzmann-factor enhanced like-charge repulsions. 

The predictions for T*{X) and p*(A) for a 1:1 electrolyte are embodied in Fig. 4. The 
asymmetry variable uj{X) = (1 — A)^/(l-|-A^) is convenient since it respects the symmetry 
A ^ 1/A. Note that cj(2) = 0.20, cu(3) = 0.40, and w(4.79) = 0.60, while the point charge 
limit (A = oo) gives u = 1. 

In Fig. 4 no hard-core terms have been included in the free energy Thus for the RPM 
{X = 1, uj = 0) one has T*^^ ~ Te ~ 0.0625 which may be compared with simulation 
estimates yielding T* ~ 0.049 (see The DH prediction for the critical density is very 
low, namely, pf^ = I/GAtt ~ 0.00497 [|; however, this increases to around p* ~ 0.03 when 
Bjerrum ion pairing is introduced |0J^. Inclusion of hard-core effects, say via a Carnahan- 
Starling form, reduces all these parameters by a few percent — and the same is expected 
to happen for the ADH-based theories. 

We note immediately from Fig. 4 that the original DH theory (i) predicts that both T*[X) 
and p*(A) rise rapidly with A. These are precisely the trends found by the MSA (using the 
energy route) [|l^ and by the MSA with Bjerrum- Ebeling-Grigo pairing |]13[. The modified 
Poisson-Boltzmann approximations of likewise predict that p* increases. (See also Fig. 
|l^.) In these approximations, however, the initial A = 1 values are well known to be 



m 



significantly higher [ T*{1) ~ 0.08, p*(l) — 0.015 to 0.03 ]; nevertheless, the proportionate 
rate of increases are roughly comparable. 

By contrast, both of the modified ADH theories predict a strong decrease in T* and p* 
when A increases from unity: see plots (ii) and (iii). Furthermore, these decreases are in 
accord with the simulations fT^JTSD (which, however, start from T*{1) ~ 0.049 and p*(l) — 



0.07) and the relative rates of fall are quite comparable. For the MF-ADH theory (ii) the 
critical temperature decreases monotonically and we find T*{X = oo) ~ 0.049 in the point 
charge limit (a++ —>■ 0), while pc exhibits a shallow minimum at A ~ 4 and then increases 
to p*(A = oo) ~ 0.006. By comparison, extrapolation of the simulations (beyond A ~ 6) 
suggests, roughly, r;(oo) ~ 0.022 and p*(oo) ~ 0.015 fll. 

In the case of the EXP"^ MF-ADH theory (iii), however, while T*{X) falls monotonically 
to about 0.053 when A — >■ oo, the critical density undergoes a shallow minimum around 
A ~ 1.8 and then rises. Insofar as the simulations seem trustworthy, and exhibit plots which 
curve downwards (i.e., are concave) vs. co'(A), it is surprising that the use of the EXP< choice 
for JF< leads to apparently inferior predictions. Indeed, on a priori theoretical grounds, the 
latter would seem superior to the MF assignment JF< = 1. We emphasize again, therefore, 
that various steps in our analysis appear most soundly based when A is not too large. 

For completeness, we report that the MF-ADH and EXP^MF^ADH theories predict 
that the critical ratio Zc = pd Pc^bT falls monotonically with A from 0.09036 at A = 1 0, 
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while the original DH theory (Sec. II. C) predicts a fall of just a few percent followed, for 
A ^ 2, by a monotonic rise. Inclusion of hard-core terms in the free energy increases these 
values (by about 7 %) but otherwise the behavior remains similar. 
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V. SUMMARY AND CONCLUSIONS 



We have extended Debye-Hiickel (DH) theory to asymmetric two-component hard-sphere 
electrolytes (i.e., "primitive models") and computed the predicted critical temperatures and 
densities. We also have derived, in Appendix A, a lower bound on the internal energy 
that extends Onsager's bound and depends only on the product of the valences 
and the unlike "collision diameter" a^_. In order to extend the original DH theory 
(based on the Poisson-Boltzmann equation) to the case of size-asymmetric ions with, say, 

a > CI++, we have identified "border zones" around ions of both species, which prove to 

be of essential importance. These zones are charge-unbalanced even at infinite temperature 
because the larger (negative) ions are geometrically excluded while the smaller (positive) 
ions may always enter: see Fig. 1. 

DH extensions which describe the border zones in a physical way (Sec. Ill) prove success- 
ful in matching trends — as determined by two independent simulation studies [n,r3| — in 
the critical temperature and density with increasing size asymmetry (see Fig. 4 and Sec. IV). 
This contrasts favorably with other theories, including several based on the mean spherical 
approximation [p!l|-p!3[|, which predict trends opposite to those revealed unequivocally by the 
simulations. 

The existence of the zones complicates the theory in an essential way; however, the usual 
DH approach can be extended straightforwardly and yields explicit approximations for the 
internal energy (and, thence, results for other thermodynamic properties). This asymmetric 
DH (or ADH) theory reproduces the limiting laws and provides the exact high-temperature 
second- virial coefficients, i?2,i and i?2,2 [see (|2.25|) - (|2.27|) ], down to the point-ion limit. 

However, in contrast to the standard DH theory for the symmetric restricted primitive 
model (RPM) with A = a = 1, the straightforward ADH theory violates the (ex- 
tended) lower bound on the internal energy in the coexistence region when A > 5. Even 
more seriously, for moderate asymmetries and moderate temperatures, the mean charge in 
a "sub-border" zone (that surrounds a -|- ion) is predicted to change sign and, at low T, 
to diverge; but, by construction of the model, such behavior is physically impossible! This 
pathology is readily traced to use of the standard DH linearization of the Boltzmann factor 
within the border zones: see Sec. II. B. Modifications of the ADH theory are thus essential 
for applications at low T. 

As shown in Sec. HI, one may restore physically sensible behavior while retaining the 
exact high-T behavior by introducing "border zone factors", J-'^(T) and T^{T), originally 
proposed in . (This also ensures that the energy bounds are no longer violated in the crit- 
ical region and beyond.) The simplest such modification amounts to a mean-field approach 
in which the Poisson-Boltzmann factors in the border zones (only) are merely replaced by 
their T = oo limits, namely, unity. For a sub-border zone (around a smaller + ion) a theo- 
retically preferable approach, dubbed EXP^, replaces the self-consistent Poisson-Boltzmann 
factor by a mean bare-interaction Boltzmann factor [or, better, linearizes about such a mean 
value: see ( ^.lUI )]. A corresponding exponential treatment of the super-border zone (around 
a negative ion) proves, however, more problematical owing to physically crucial charge- 
saturation effects that are hard to elucidate in a precise way. One may expect that both the 
mean-field and EXP*~ modifications of linear ADH theory are most reliable at small degrees 
of asymmetry. 
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To explore the consequences of the simplest modified mean-field and EXP^ ADH theories 
(just sketched) we have computed the predicted critical parameters as a function of the size- 
asymmetry A for 1:1 electrolytes [with additive interactions, a+_ = |(a++ + a )]: see Fig. 

4. As A increases from unity, the predicted T*{X) and p*(A) fall systematically within both 
of these modified ADH theories. [See ( |1.4| ) and ( |1.5| ) for definitions of the reduced units.] 
These decreases accord well with recent simulations |ll^ , p^ . On the other hand, an original 
proposal by Debye and Hiickel in 1923, that completely ignores the border zones (see Sec. 
III.C), predicts diametrically opposite trends. Furthermore, current, more sophisticated 



theories |]T^,|T3|] make similar predictions of increasing T* and p* (in addition to yielding 
excessively large values of T* for the RPM ^j^). 

We conclude that DH-based theories seem to extract the basic physics in a quantitatively 
more reliable way even for size-asymmetric systems, than do potentially better, but 

physically less transparent approaches like the MSA. It is still necessary, however [|l^], to 



incorporate Bjerrum ion-pairing and dipole-ion solvation [^ into the modified ADH theories 
expounded here. (This will also increase the predicted critical densities to better match the 
anticipated values.) It is not obvious that the correct trends with asymmetry (accepting the 



validity of the simulations |T^JT^) will survive these extensions: it seems likely, nonetheless, 
that the proper dependence on asymmetry will be preserved (as suggested, e.g., by comparing 
the results of |12[ and [|13l). 

From a broader perspective, it remains frustrating that more powerful and definitive 
theoretical techniques have not yet been devised to aid in our understanding of such a 
fundamental and significant model for condensed matter. 



ACKNOWLEDGMENTS 

We appreciate the interest, encouragement, and comments of Professors Benjamin P. 
Vollmayr-Lee, Elhot Lieb, and A. Z. Panagiotopoulos. We also are grateful to Dr. S. Banerjee 
for pointing out an oversight in the appendix. Professor Katharina Vollmayr-Lee generously 
translated sections of [|l| related to the present work, and Youngchan Kim kindly commented 
on a draft typescript. The bulk of the research reported here was supported by the National 
Science Foundation (under Grant Nos. CHE 96-14495 and CHE 99-81772). D.M.Z. is grateful 
for current support from Professor Thomas B. Woolf and the National Institutes of Health. 
S.B. gratefully acknowledges the current support of Professor Terry Gaasterland. 



18 



APPENDIX A: ENERGY BOUNDS FOR THE TWO-SPECIES PRIMITIVE 

MODEL 



We here generalize the Onsager bound for the RPM to the general two-component 
hard-sphere model defined in Sec. II and also develop an (apparently new) lower bound 
which improves on the simple Onsager result. 

The quantity to be bounded is the interaction energy of the ions, [/^°*. Defining the 
interionic distances, r,- 



ij , and interaction potentials, 



(Al) 



where is the charge on the ith ion, one has 



i<j 



(A2) 



By adopting the 1 / r jj form in (|A1|) , we implicitly assume that the charge densities between 
any two ions are spherically symmetric and never overlap. To avoid irrelevant singulari- 
ties, we further assume that the charge density is everywhere finite, i.e., has been suitable 
"smeared" or distributed. 

consider the total electrostatic energy, [7*°*, which, besides the 



Following Onsager [18 



interaction energy, also includes the self energies, , i.e., the energies required to assemble 
the individual ions from an infinitely dispersed charge state of zero energy. As is well known 
(see, e.g., |^), the total energy may be expressed in terms of the electric field, E, as a 
positive definite quantity, namely, 



U 



tot 



i<3 



self 



— /d^r |E|^ > 0. 

Stt 



(A3) 



Note that the smeared charge distributions assumed guarantee the finiteness of [/*°*. 



1. Generalization of the Onsager Bound 

The Onsager bound for the RPM and its direct generalization follow easily from ( |A3D . 
Rearrangement leads to 

i<j i 

This inequality holds for arbitrary charge distributions obeying the restrictions stated above; 
different distributions, of course, lead to different values of uf^^. To obtain the strongest 
bound, we must minimize the magnitudes of the self energies. This is accomplished by 
dispersing the ionic charge as much as possible, namely, by placing it in a thin shell on the 
surface of the largest permissible sphere (of diameter a™'^^). Considering a vanishingly thin 
shell, the minimal self energy is thus 

min {uf^ = g^V^^r"- (A5) 
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For the symmetric RPM, the equahty of charges, g+ = — g_ = g, and of sizes, af^^^ = = 

a, immediately leads to the Onsager bound q^/Da. In the case of non-additivity, however, 
one must avoid overlapping charge distributions from distinct ions, so that 

^max ^ {a+^, an} . (A6) 

We can now formulate the explicit Onsager bound. First, defining z = \z+/z^\, note that 
the total numbers of particles of each species are 

= N/{1 + z) and A^. = zNjiX + z) . (A7) 

Recalling the definition of the reduced energy ( |1.6| ) and using the constraint of overall charge 
neutrality, one finds the bound, namely, 

m(p, T) > Mons = — • (A8) 

L ~\- Z 

Note that for the size-symmetric case, where the charge asymmetry does 

not affect the bound in these reduced units, so that u(p, T) > wons = — 1- If we separate the 
three basic size asymmetry cases, the result ([A8|) translates into 



«Ons = - W[(l + + A\ + + + m , for a++, a__ < a 

«Ons = - + + ^)] + 1/(1 + A} , for a++ < a < a__ (A9) 

Mons = ~1 ) for a < a+_|_, a . 

Note that the strongest bound is -1, which obtains for Case III, when the like-diameters 
both exceed a, which matches the RPM result. On the other hand, wons is weaker in Cases I 
or II, since the bound can then diverge to — oo when b„ ^ —\ (or, equivalently, as Oo-o- 0). 
Physically, shrinking the like-diameters (but keeping a+_ positive) should not decrease the 
energy: thus a stronger bound is desirable. 



2. Improved Bound 

To do better we compare the two-species primitive model with another model whose 
energy can be bounded by an Onsager construction, specifically, an interpenetrating, two- 
species shell model, consisting of a charge-neutral mixture of uniformly surface-charged 
spheres of total charges = 2:4.^0 or z_qi^ and equal diameters a, but with no hard-core 
constraint. The shell diameter a will be identified with a+_ for the primitive models. The 
interaction potential is described further below. 

Since the configuration space of the primitive model — in terms of ion-center locations 
— is a subset of the space of the shell model, any ground state configuration of the primitive 
model of energy, say f/™*, is present in the shell model with energy, say f/™*, which cannot 
be lower than the shell model ground state, say f/"g. Thus, if we can establish 

f/^°* > t/f* (AlO) 
and also show that f/^g is bounded below, we will obtain a lower bound on f/™*, as desired. 
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We may construct the interaction potential of the interpenetrating shell model, say iptj 
(which does not behave as l/r^j at all separations), as the difference between the total energy 
of a two-particle system and the sum of the two self energies. Thus if /)fc(r;rfc) represents 
the charge distribution of a shell ion k centered at r^, and 

p(r) = p,{r; r,) + pj{r; rj) (All) 

is the total charge density for two ions, the interaction potential follows from 

where E(r; r^, r^) and 0(r; r^, r^) are the total field and electrostatic potential. If we now 
define (/)fc(r; r^) to be the electrostatic potential resulting from an isolated shell ion k at r^, 
we have 

2 

i y'dVpfc(r;rfc)0,(r;rfc) = . (A13) 



Then, using the linearity of the charge density and potential, the relation ( [A12| ) may be 
simplified to yield 

^ij = <P{rij) = 1 Jd^r [pi{r; ri)(f)j{r; rj) + pj{r; rj)0»(r; r^)] . (A14) 

To compare the energy of an arbitrary primitive-ion configuration with that of the corre- 
sponding shell- ion configuration in order to establish (|A10|) , we observe first that because of 
the pairwise additivity of the interaction energy ( [A2[ ), one need analyze only two shell ions 
of the same charge which overlap. To see this, note that all (+, — ) ion pairs in a primitive 
model are separated by distance not less than = a, which is the same diameter as the 
shell ions. Thus, oppositely charged shell ions that correspond to a primitive ion configu- 
ration never overlap. The only differences arise when, in the primitive ion system, one or 
both of the like diameters, a„„, is smaller than a. This will allow overlapping shells in the 

corresponding shell configurations. If a^^ and a exceed a, the energy of corresponding 

primitive and shell configurations will always be identical. 

Consider then, for concreteness, two positive overlapping shells, separated by a distance 
r++, with a++ < r++ < a; see Fig. 5. We want to show that the interaction energy of 
this pair, (^(r++), is less than that for the the corresponding non-overlapping primitive ions 
which is simply q'^/Dr^^. 

By symmetry the two terms in (|A14D are now identical. Thus, consider the charge 



distribution of the right-hand shell in Fig. 5 in the potential of that on the left; the right-hand 
distribution divides naturally into the two parts shown in the figure: a part exterior to the 
left-hand shell (bold) and an interior part (dashed). If r; denotes the position of the left-hand 
ion, the resulting potential at an exterior points, r>, is simply q+/\ri — r>| (because these 
points "see" a spherically symmetric left-hand charge distribution). Conversely, interior 
points, such as r<, experience only the constant electrostatic potential g+/(a/2). This is 
clearly less than the potential they would experience were all the left-hand charge distributed 
on the smaller primitive ion sphere with the same center, r^. Consequently, overlapping like- 
charged shell ions have a smaller interaction potential than the corresponding primitive ions 
with the same centers. This establishes ( |A1CI|) . 
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To obtain a lower bound for the shell model itself, recall ( [X^) and ( |A11| ) and bound the 
total energy of any shell configuration as 

= i- Jd^r |E|2 = 1 yd\p(r)0(r) > , (A15) 

where the total shell charge density, p(r) can be expressed in the form ( [Allj ) but with a sum 
extending over all the shell ions. The total electrostatic potential, 0(r), can be decomposed 
similarly, yielding 

\jd'rp{v)(t){v) = [pi(r) + ■ ■ ■ + p^iv)] [01 (r) + ■ ■ ■ + 0^(r)] . (A16) 

Finally, by combining the previously defined shell self energies, ( |A13|) , and interaction po- 
tentials in ( |A14|) , the inequality ( |A15| ) may be re-arranged to give 

^ = lv.3..>^>-l + ^\ = (Km 

N N "^.^'^ - N - N \ Da Da Da ' ^ ^ 

The bound on the primitive model is now completed by combining this with (|A10|) . In 
terms of the reduced energy per particle, (p. .61), the result may be written 



u{p,T)>-l. (A18) 

Note that the positive definite collision diameter, a^ = a, and the valences Za do not appear 

explicitly here since they enter into the definition ( |1.6|) . Since a++ and a are also absent, 

the bound remains valid for point ions {acra 0). 
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Figure Captions 



1. Illustration of a border zone. The grey sub-border zone may be occupied only by the 
centers of + ions since a+_)_ < a+_. 

2. Effects of size asymmetry on high- and low-temperature energy isochores, according to 
the linear ADH theory for a 1:1 primitive model. For each density, the RPM isochores 

(with A = a = 1) are shown as solid curves. The associated, successively lower 

plots correspond to A = 2 (dashed), 4 (dot-dashed), and 6 (dotted, for p* = 0.1 only). The 
density-dependent, T* —* oo limiting values are exact. (The near-agreement between two of 
the curves at T* = oo is a coincidence.) 

3. Effects of size asymmetry on the low-temperature energy isochores, according to the 
simplest "mean-field" modification of ADH theory for a 1:1 primitive model at densities 
p* = 0.01, 0.03, and 0.1 and size asymmetries A = 1 (the RPM) solid curves, A = 2, dashed, 
and A = 4, dot-dashed curves. 



4. Critical temperature and density predictions for a 1:1 electrolyte with additive hard- 
sphere interactions as a function of the size asymmetry variable u){X) = (1 — A)^/(l -|- A^), 

which increases monotonically with A = a The reduced parameters p* and T* 

are defined via (|1.4|) and (|1.5| ): (i) represents the 1923 DH theory, while modifications of 
asymmetric DH theory are embodied in (ii) with mean-field factors for both J-'^ and 
and (iii) with a Boltzmann-factor (EXP^) border zone factor for JF< and a mean- field factor 
J^^ . The circles denote the T* simulation estimates of Romero- Enrique et al. WM. 



5. Two overlapping shell ions, with corresponding hard-core primitive ions (dotted), which 
are smaller but concentric. The shell diameter, a, is identified with a+_ of the primitive 
models. 
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FIGURES 




FIG. 1. Zuckerman, Fisher, and Bekiranov 



25 




FIG. 2. Zuckerman, Fisher, and Bekiranov 
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FIG. 3. Zuckerman, Fisher, and Bekiranov 
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FIG. 4. Zuckerman, Fisher, and Bekiranov 
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